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Circadian rhythms in mammals are controlled by the neurons located in the 
suprachiasmatic nucleus of the hypothalamus. In physiological conditions, the sys- 
tem of neurons is very efficiently entrained by the 24-hour light-dark cycle. Most of 
the studies carried out so far emphasize the crucial role of the periodicity imposed 
by the light dark cycle in neuronal synchronization. Nevertheless, heterogeneity as 
a natural and permanent ingredient of these cellular interactions is seemingly to 
play a major role in these biochemical processes. In this paper we use a model that 
considers the neurons of the suprachiasmatic nucleus as chemically-coupled modi- 
fied Goodwin oscillators, and introduce non-negligible heterogeneity in the periods 
of all neurons in the form of quenched noise. The system response to the light-dark 
cycle periodicity is studied as a function of the intcrncuronal coupling strength, 
external forcing amplitude and neuronal heterogeneity. Our results indicate that 
the right amount of heterogeneity helps the extended system to respond globally in 
a more coherent way to the external forcing. Our proposed mechanism for neuronal 
synchronization under external periodic forcing is based on heterogeneity-induced 
oscillators death, damped oscillators being more cntrainablc by the external forcing 
than the self-oscillating neurons with different periods. 

Keywords: Circadian oscillations; quenched noise; noise-induced oscillators 
death; modified Goodwin model; noise-induced synchronization. 



1. Introduction 

Circadian rhythms are light-dark dependent cycles of roughly 24 hours present 
in the biochemical and physiological processes of many living entities (Rcppert & 
Weaver 2002). In mammals the main mediator between the light-dark periodicity 
and the biological rhythms is formed by two interconnected suprachiasmatic nuclei 
(SCN), located in the hypothalamus. These nuclei form the so called "circadian 
pacemaker" and contain about 10.000 neurons each (Reppert & Weaver 2002; Moore 
et al. 2002). 

The main property of the SCN is that their activity displays self-sustained oscil- 
lations in synchrony with the external forcing imposed by the light-dark cycle. The 
exact mechanism leading to this behavior has been the subject of intense research. 
It has been shown that, when taken individually, neurons produce oscillations with 
a constant period ranging from 20 to 28 hours (Honma et al. 2004; Welsh et al. 
1995). The oscillatory behavior originates in a regulatory circuit with a negative 
feedback loop. The relevant question is how this individual oscillatory behavior 
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translates into common, global, oscillations of the SCN activity synchronized with 
the external light stimulus. 

It has been shown that the origin of the oscillatory activity of the circadian 
pacemaker at the global level resides on the interaction between the SCN neurons. 
Coupling between cells in the SCN is achieved partly by neurotransmitters (Honma 
et at 2004; Hastings & Hcrzog 2004) and it is by means of those neurotransmitters 
that external forcing by light influences the neuronal synchronization. For example 
the vasoactive intestinal polypeptide (VIP) has been shown to be necessary in medi- 
ating both the periodicity and the internal synchrony of mammalian clock neurons 
(Shen et al. 2000; Aton et al. 2005; Maywood et al. 2006). Therefore, a model of 
coupled and forced neurons appears quite naturally as responsible for the circadian 
rhythms. Along these lines, an interesting mechanism has been put forward recently 
by Gonze et al. (2005) and by Bernard et al. (2007). They proposed that synchro- 
nization to the external forcing is facilitated by the fact that interneuronal coupling 
transforms SCN into damped oscillators which can then be easily entrained. 

In this paper we show that the presence of some level of heterogeneity or dis- 
persion in the intrinsic periods of the oscillators (Schaap et al. 2003; Herzog et al. 
2004) can improve the response of the coupled neuronal system to the external 
light-dark forcing. The proposed mechanism for the improvement of the neuronal 
synchronization under external periodic forcing bears some similarities with the one 
proposed in (Gonze et al. 2005; Bernard et al. 2007) in the sense that the oscillators 
are brought to a regime of oscillator death (Ermentrout 1990; Mirollo & Strogatz 
1990), but in our case this regime is induced by the presence of heterogeneity. Once 
this regime has been reached, the damped oscillators are more entrainable by the 
external forcing than the self-oscillating neurons with different periods, or the syn- 
chronized oscillatory state which appears in the strong coupling regime but with a 
period larger than the individual neuronal periods. 

To be more specific, we will assume that the periods of the individual neurons 
are random variables drawn from a normal distribution. We will then analyze the 
global response of the system to the light-dark cycle periodicity as a function of the 
interneuronal coupling strength, external forcing amplitude and neuronal hetero- 
geneity. We show that the presence of the right amount of dispersion in the periods 
of the neurons can indeed enhance the synchronization to the external forcing. 

Period dispersion arises as a consequence of the cellular heterogeneity at the 
biochemical level, which is an experimentally well observed fact (Aton & Hcrzog 
2005; Honma et al. 2004). It can act in cither physiological or pathological condi- 
tions. An example of the latter is the diversification of antigenic baggage present 
in tumor cells that makes them more difficult to be recognized and captured by 
the defense mechanisms and therefore more prone to migrate and develop metasta- 
sis (Gonzalez-Garcfa et al. 2002). Our results show that some level of disorder can 
be of help when synchronizing neuronal activity to the external forcing. Although 
counterintuitive, it has been unambiguously shown that the addition of various 
forms of disorder can improve the order in the output of a large variety of non- 
linear systems. For example, the mechanism of stochastic resonance (Gammaitoni 
et al. 1998; Hanggi & Marchcsoni 2009) shows that the response of a bistable sys- 
tem to a weak signal can be optimally amplified by the presence of an intermediate 
level of dynamical noise. Stochastic resonance is not a rare phenomenon; it has been 
repeatedly shown to be relevant in physical and biological systems described by non- 
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linear dynamical equations (Gammaitoni et al. 1998; Hanggi & Marchcsoni 2009). 
In large systems with many coupled elements, noise is responsible for a large variety 
of ordering effects, such as pattern formation, phase transitions, phase separation, 
spatiotemporal stochastic resonance, noise-sustained structures, doubly stochastic 
resonance, amongst many others (Garcfa-Ojalvo & Sancho 1999). All these exam- 
ples have in common that some sort of order at the macroscopic level appears only 
in the presence of the right amount of noise or disorder at the microscopic level. 
Furthermore, it has been proven that noise may play a constructive role in nonlin- 
ear systems, by enhancing coherent (periodic) behavior near bifurcations and phase 
transitions (Neiman et al. 1997; Pikovsky & Kurths 1997). In this paper we intro- 
duce non- negligible random heterogeneity into the periods of all neurons, so-called 
quenched noise. Numerical simulations suggest (data not shown) that the results arc 
valid as well when the quenched noise is introduced into the model parameters. A 
different approach is the consideration of intracellular stochastic variability due to 
low molecule numbers (Forger & Peskin 2005) or both variability and heterogeneity 

Close to our work is the study by Ueda et al. (2002), where the effect of fluc- 
tuations in neuron parameter values is assessed and it is shown that the coupled 
system is relatively robust to noise. Previous theoretical studies have addressed the 
effect of noise on genetic oscillators (Thattai & van Oudenaarden 2001; Steuer et al. 
2003; Becskei et al. 2005), and some have proposed an ordering influence of noise 
on circadian clocks at the single cell level in cases where neither light intensity nor 
coupling strength by themselves can synchronize the system. Collective phenomena 
induced by heterogeneity in autonomous, non-forced systems, has also been dis- 
cussed in the literature. For example (de Vries & Sherman 2001) and (Cartwright 
2000) have shown that collective bursting or firing can appear in excitable systems 
and a general theory of the role of heterogeneity in those systems has been devel- 
oped by (Tessone et al. 2007). In this paper, we refer to the collective response 
in systems of non-linear oscillators subjected to the action of an external forcing 
representing the day-light cycle. 

The paper is organized as follows. In section 2 we will describe in detail the 
model of circadian oscillators and the methods we use. It is a coupled extension of 
the original Goodwin oscillator (Goodwin 1965) as developed by Gonze et al. (2005). 
In section 3 we analyze the system response to the periodic external forcing, as a 
function of the external forcing amplitude, coupling strength and neuronal diversity 
or heterogeneity. By simulating numerically the governing differential equations we 
identify the range of these parameters for which the extended system oscillates in 
synchrony and entrained to the external light period. Section a describes the mech- 
anism through which the neuronal heterogeneity favors the synchronization with 
the external forcing and analyzes the combined influence of the coupling strength, 
neuronal heterogeneity and light amplitude on the stability of the linearized sys- 
tem of coupled oscillators. We show that a mean variable in this model exhibits 
a transition from a rhythmic to an arrhythmic dynamics (the so-called oscillator 
death (Ermentrout 1990; Mirollo & Strogatz 1990)). Concluding remarks arc found 
in section 4. 
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2. Model and methods 

(a) The circadian pacemaker 

As stated in the introduction, our aim is to consider the role that the hetero- 
geneity in the population of neurons plays in the global response of the SCN to 
an external oscillating stimulus. To this end, we consider an ensemble of coupled 
neurons subject to a periodic forcing. Each of the neurons, when uncoupled from 
the others and from the external stimulus, acts as an oscillator with an intrinsic 
period. Heterogeneity is considered insofar the individual periods are not identi- 
cal, but show some degree of dispersion around a mean value. For each one of the 
neurons in the SCN we use a four-variable model proposed by Gonze et al. (Gonzc 
et al. 2005), which is based originally on the Goodwin oscillator (Goodwin 1965), 
to describe circadian oscillations in single cells. The variables of the model are as 
follows: The clock gene mRNA (X) produces a clock protein (Y), which activates a 
transcriptional inhibitor (Z) and this in turn inhibits the transcription of the clock 
gene, closing a negative feedback loop. The mRNA X also excites the production 
of neurotransmitter V , which in the coupled system will be then the responsible 
of an additional positive feedback loop. In order to overcome the high Hill coeffi- 
cients required for self-oscillations, Gonze et al. replaced the linear degradation by 
nonlinear Michaclis-Mcntcn terms. This leads to the system of equations: 

dX K{ X , n . 

= V\ — — 7-^2 , (2.1) 

dt Kf + Z 4 A 2 +A' V ' 

dY Y 

IT = ksX-vtjrr-v, (2.2) 
dt K4 + Y 

§ = k 5 Y-, 6 -^—, (2.3) 
dt K 6 + Z 

at K s + V 

which, depending on parameters, might produce oscillations in a stable limit cycle. 
Using the values v\ = 0.7 nM/h, v 2 = Va = v§ = 0.35 nM/h, v% = 1 nM/h, 
A'i = K 2 = Ki = K 6 = A' 8 = 1 nM, fc 3 = k 5 = 0.7/h, k 7 = 0.35/h, the period of 
the limit cycle oscillations is T = 23.5 h. 

For the complete model, we take N neuronal oscillators, each one of them de- 
scribed by four variables (Xi, Y, Zi,Vi), i = 1, . . . , N , satisfying the above evolution 
equations. Heterogeneity in the intrinsic periods is introduced by multiplying the 
left-hand-side of each one of the equations (2.1-2.4) by a scale factor Tj. Hence, 
the intrinsic period Tj of the isolated neuron i is t{T . The numbers Tj are indepen- 
dently taken from a normal random distribution of mean 1 and standard deviation 
<7. Since the periods must be positive, in the numerical simulations we have ex- 
plicitly checked that, for the values of a considered later, r, never takes a negative 
value, which would be unacceptable. The standard deviation a will be taken as a 
measure of the diversity. A value of a = 0.1 for example corresponds to a standard 
deviation of 10% in the individual periods of the uncoupled neurons, close to the 
observed variation of periods between 20 and 28 hours. 

Two additional factors influence the dynamics of single cell oscillations: forcing 
by light and intercellular coupling. Both are assumed to act independently from the 
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negative feedback loop and arc added as independent terms in the transcription rate 
of X (Gonzc et al. 2005). Light is incorporated through a periodic time-dependent 
function L(t), which can be realized in various forms. In the majority of the pre- 
sented results we have used a sinusoidal signal, L(t) = (1 + sinw£). In some 
cases, for comparison and to simulate different day lengths, we have used a square 

wave L(t) = < ^ , mod 24ft) < tu g ht ^ j n ways the signal oscillates be- 

w \ 0, otherwise 6 
tween the values L(t) = and L(t) = Lq with a period 2tt/u> = 24ft. 

Coupling between the neurons is assumed to depend on the concentration F of 
the synchronizing factor (the neurotransmitter) in the extracellular medium, which 
builds-up by contributions from all neurons. Under a fast transmission hypothesis, 
the extracellular concentration is assumed to equilibrate to the average, mean-field, 
cellular neurotransmitter concentration, F = The resulting model is: 

Ts = v\ — -, — - — 7 — ^2 Vv c h Lit) (2.5) 

1 dt l Kf + Zf K 2 +Xi K C + KF w 1 ' 

dY Y- 
r~ = k^-^^-^r, (2.6) 

= kM-Vs-^-, (2.7) 
dt K 6 + Zi 

dVi T, V V i tO Q\ 

dt A 8 + Vi 

N 



F 



N 

i=l 



with v c = 0.4 nM/h, K c = 1 nM. 

There is experimental evidence supporting the assumption of a chemical (rather 
than electrical) mechanism of inter-cell communication among SCN neurons as a 
synchronization factor and, in fact, mechanisms other than neurotransmitters or 
electrical coupling for the SCN communication have been suggested (e.g. by Pol & 
Dudck (1993)). Furthermore, more realistic modeling which takes into account all 
variables known to participate of the negative feedback loop has been introduced. 
These models may include up to 10 variables and corresponding equations for each 
single cell (Bernard et al. 2007). 

It seems, however, that in order to get understanding of the SCN dynamics, a 
sufficient tool is the 4 variable model described above. In fact, the synchronization of 
damped oscillators is independent from the particular intracellular model used and 
as discussed by (Bernard et al. 2007), this system, the model developed by (Leloup 
& Goldbcter 2003), and other simple negative feedback oscillators have similar 
synchronization properties. In this paper we have decided to use the simpler 4- 
variablc model although most of our results are also valid in the more complex 
10- variable model. 

A model close to (2.5-2.9) has been used by Ullner et al. (2009), where the 
authors investigate how the interplay between fluctuations of constant light and in- 
tercellular coupling affects the dynamics of the collective rhythm in a large ensemble 
of non-identical, globally coupled oscillators. In their case, however, an inverse de- 
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pendence of the cell-cell coupling strength on the light intensity was implemented, 
in such a way that the larger the light intensity the weaker the coupling. 

(b) Measures of synchrony and entrainment 

Due both to the effect of coupling and of forcing, the neurons might synchronize 
their oscillations. There are several possible measures of how good this synchroniza- 
tion is. In this paper, the interneuronal synchronization will be quantified by the 
parameter of synchrony p, defined as 



1 - 



(2.10) 



where (...) denotes a time average in the long-time asymptotic state. The pa- 
rameter p varies between a value close to (no synchronization) and 1 (perfect 
synchronization, with all neurons in phase, = Vj(t),Vi,j). It is important to 
note that even if the neurons synchronize perfectly their oscillations, the period 
of those oscillations does not necessarily coincide with the mean period T of the 
individual oscillators or with the period 2-k/uj of the external forcing. In fact, in 
the unforced (no light) case, the period of the common oscillations (for the set of 
parameters given before and a dispersion of a = 0.05 and coupling K = 0.5) is 
approximately equal to 26.5 h whereas the period of the forcing is 2tt/uj = 24 h 
and the mean period of the individual uncoupled oscillators is T = 23.5 h (Gonzc 
et al. 2005). 

Besides the previous measure of synchronization amongst the oscillators, we are 
also concerned about the quality of the global response of the neuronal ensemble to 
the external forcing L(t). A suitable measure of this response can be defined using 
the average gene concentration, 



1 N 

x(*) = ]v£*i(*) 



(2.11) 



and computing the so-called spectral amplification factor R (Gammaitoni et al. 
1998), 

- lult X(t))\ 2 . (2.12) 



R =j2 



R is nothing but the normalized amplitude of the Fourier component at the forcing 
frequency u> of the time series X(t). We will show that, under some circumstances, 
the response R will increase with the intrinsic diversity a and that the period of 
the oscillations at the global level coincides with that of the external forcing, these 
being the main results of this paper. 



3. Results 

The synchronization properties of the set of circadian oscillators is influenced by 
the amplitude of the external forcing Lq , the coupling strength K and the diversity 
in the individual periods a. The role of the first two has been studied in (Bernard 
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Figure 1. Synchrony order parameter p (see Eq. (2.10)). Values are coded in colour levels, 
and displayed as a function of Lo and a for several values of K. Data from numerical 
simulations of N = 1000 neurons with dynamics ruled by Eqs. (2.5-2.9). Synchrony among 
the neurons (yellow region) is favored by strong or very weak light intensity Lq, low 
diversity a and large coupling K. The thick black line is the linear stability limit discussed 
in section a (see also Fig. 7). 

et al. 2007; Becker- Wcimann et al. 2004; Gonze et al. 2005). In this section we focus 
on the heterogeneity of neuronal periods and analyze the combined influence of Lq, 
K and a on the different parameters quantifying interneuronal synchronization and 
response to the forcing. 

Fig. 1 shows colour plots of the parameter of synchrony p as a function of the 
diversity a and the light intensity Lq, for different values of the coupling strength 
K. High values of the light intensity Lq favor interneuronal synchrony. Also in 
agreement with its intuitive disordering role, high neuronal diversity leads to a low 
synchrony parameter p in several parts of the diagrams. However, there is a region 
of values of Lq € [0, L max ] for which there is a non-monotonous dependence of the 
synchrony order parameter with respect to the diversity. This can be seen more 
clearly in panel (a) of Fig. 2 where we plot p as a function of diversity a for fixed 
values of K = 0.6 and Lq = 0.005. p first decreases by increasing a within the 
interval < a < 0.05, but then it develops a maximum. The range of values of 
Lq for which this non-monotonous behavior is observed depends on the coupling 
constant K: the larger K, the larger L max . 

As stated before, the fact that neurons synchronize amongst themselves does 
not mean that they synchronize to the forcing by light. To study this point, we have 
computed the individual periods TJ, i = 1, . . . , N, of the oscillators in the ensemble. 
In those cases in which the concentrations do not oscillate with exact periodicity, 
we still define the period as the average time between maxima of the dynamical 
variables. In Fig. 4 we plot the mean value T = as a function of a and 
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Figure 2. Main parameters used for characterizing the synchronization of circadian oscil- 
lators as a function of the variance a. (a) the synchrony parameter p; (b) the mean T of 
the individual periods Tt; (c) the response order parameter R; (d) the maximum real part 
of the eigenvalues of the linearized system. 

Lq for different values of K. As the dispersion in Tj is small, it turns out that T is 
close to the period of the average variable X(t). 

Although, by construction, individual neurons have periods that fluctuate around 
T = 23.5 h, it turns out that the period of the resulting synchronized oscillations 
that occur in the unforced but coupled {Lq = 0, K > 0) case, increases with in- 
creasing coupling K. For example, T ps 30 h for K = 0.6, mostly independent of the 
value of a. As the forcing sets in, at low values of the coupling strength, the mean 
period is now T = 24 h for all values of Lq and a. As the coupling between neurons 
increases, larger values of Lq and/or a arc needed in order for the mean period to 
coincide with that of the external forcing. An important feature that emerges from 
these plots is that for low light intensity it is possible to achieve a mean period 
of 24 h by increasing the neuronal diversity. For example, in the areas at the left 
of the different panels of Fig. 4, or in panel (b) of Fig. 2 corresponding to the case 
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Square wave light (K=0.6, L =0.005) 
8L:16D 12L:12D 16D:8L 




Figure 3. Characterization of synchrony for a light signal of square wave form. The syn- 
chrony parameter p, the mean T of the individual periods T; and the response order pa- 
rameter R (from top to bottom) measured for square wave stimuli of various day lengths 
(8/1, 12h and 16h, from left to right). 

K = 0.6, while identical neurons have periods of w 30 h, increasing a induces an 
adjustment of the period to 24 h. The transition between T — 24 h and T ^ 24 h 
is rather sharp, specially for large K . This is a clear manifestation that diversity 
indeed is able to improve the response to the external forcing. The same conclusion 
about the constructive role of diversity can be reached by looking at the measure 
of response R (see Figs. 5 and 2(c)). These figures show that there is a region in 
parameter space in which the system response to the periodic light forcing displays 
a maximum value as a function of diversity a. This indicates that it is possible 
to improve neuronal synchronization to the daily-varying light input by taking a 
close to an optimal value. Too small or too large diversity will not yield an optimal 
response. This is a clear manifestation that diversity indeed is able to improve the 
response to the external forcing. 

A complementary perspective on this constructive role of diversity is attained 
looking at spectral amplification factor, R, from Eq. (2.12). This is a normalized 
measure of the amplitude of the oscillation of the neuronal system at the frequency 
of the daily forcing. Figures 5 and 2(c) show that there is a region in parameter space 
in which the system response to the periodic light forcing displays a maximum value 
as a function of diversity a. In fact this maximum is very large as compared with 
the R value at zero diversity, so that one can say that one of the most noticeable 
effects of a non-vanishing neuronal diversity is to give the system the capacity to 
respond efficiently to the 24h forcing in situations of small or no response at this 
frequency in the absence of diversity (the non-diverse neuronal ensemble could be 
oscillating at a different frequency, as revealed by high values of p). In summary, 
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Figure 4. Colour level plots of the mean of the individual periods T of 1000 neurons 
under forcing at 24h cycle for the system of Eqs. (2.5-2.9). High light intensity Lo and 
high diversity a assures entrainment of oscillations to external frequency (blue region). 
Increasing coupling enlarges the region (yellow) of oscillations at a period larger than that 
of the driving force. 



it is possible to largely improve neuronal synchronization to the daily-varying light 
input by taking a close to an optimal value. Too small or too large diversity will 
not yield an optimal response at this frequency, although the response is generally 
larger than for zero diversity. 

An external signal of square wave form and with different day lengths lead to 
similar results. As can be seen in figure 3 the response R to the external signal 
passes through a maximum at a intermediate value of diversity. The mean period 
and the synchrony parameter behave as in the case with a pure sinusoidal as the 
driving force. Furthermore, the qualitative result is independent of the chosen day 
length. 



(a) Diversity-induced oscillator death 

Why does an increase in the diversity of the oscillators lead to an improved 
response to the external forcing? We argue that the main effect of the increase of 
the diversity is to take the oscillators into a regime of oscillator death (Ermentrout 
1990; Mirollo & Strogatz 1990) in which they can be easily entrained by the vary- 
ing part of the forcing. To understand this mechanism we first split the forcing into 

a constant (the mean) and a time varying part: L(t) = — H — — sin(<Jt). Taking 

only the constant part, L(t) = Figs. 6(a)— (c) show that the oscillators go from 

self-sustained oscillations to oscillator death, i.e. the amplitude of the self-sustained 
oscillations decreases, as a increases. Once oscillators are damped, they would re- 
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Figure 5. Colour plots of the spectral amplification factor R as defined in Eq. (2.12) in 
logarithmic scale. Too high light intensity Lo and too much diversity a lower the response 
of X to the external frequency. This also happens for low light and low diversity (where 
the neurons oscillate at a frequency 7^ 24 h, see Fig. 4). Between both limits R passes a 
local maximum. 



spond quasi-linearly to periodic forcing, at least if this forcing is not too large, and 
linear oscillators always become synchronized to the external forcing, independently 
of their internal frequency. This is consistent with what is seen in figures 6(d)-(f), 
where the neurons in the case of low heterogeneity oscillate synchronously with each 
other, but their common period is larger than the one of the light forcing. Only 
when diversity brings the neurons to oscillator death can all of them be entrained 
to the period of the forcing signal. The mechanism is related to the one discussed 
by Gonze et al. (2005) and Bernard et al. (2007), but here we stress that neuron 
heterogeneity, as opposed to internal neuron parameters and couplings, is enough 
to damp the collective neuron oscillations and bring the system to a non-oscillating 
state where it can be more easily entrained. It is interesting to note that it has 
been shown experimentally for fruitflies that only a subset of the pacemaker neu- 
rons sustain cyclic gene expression after changing the laboratory light conditions 
to constant darkness, whereas the oscillations of the other pacemaker neurons are 
damped out (Veleri et al. 2003). Although this does not reveal the mechanism by 
which the oscillations die out it suggests that some of the circadian oscillators do 
indeed work in the damped regime, at least in Drosophila. 

An alternative way of checking this mechanism based on diversity-induced os- 
cillator death is by analyzing the stability of the steady state of the system of Eqs. 

(2.5-2.9) when considering a constant forcing L(t) = ——. The numerical calcula- 
tion of the fixed point of the dynamics is greatly simplified by the fact that the 
concentrations of the biochemical variables are the same for each one of the N neu- 
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Figure 6. Figures (a), (b) and (c) represent the time-dependent amplitude of the Vi variable 
for a few selected neurons in the presence of constant light and increasing a, while Figures 
(d), (e) and (f) represent the amplitude of the same neurons with sinusoidal light and 
increasing a. The thin line on the bottom of the graphs is the external light signal. K — 0.6 



rons irrespectively of their specific value of t%. The system (2.5-2.9) is linearized 
around this steady state and the eigenvalues of the stability matrix computed for 
several realizations of diversity parameters r,;. In each case, the positive or nega- 
tive character of the real part of the eigenvalue with the largest real part indicates 
the instability or stability, respectively, of the fixed point solution. In Fig. 7 we 
plot the mean of that maximum real part of the eigenvalues averaged over various 
realizations of the time scales Tj, for N = 200 coupled neurons, as a function of 
Lq and (j, and different values of the coupling K (see also panel (d) in Fig. 2). In 
every diagram we can see that low diversity or low forcing yield an unstable steady 
state (yellow region). This is where self-sustained oscillations are observed. A thick 
black line in the contour plots indicates a zero real part. The relevance of this line 
separating positive from negative maximum average eigenvalues is more apparent 
when we note that it also delimits regions of interest in Figs.l, 4 and 5. 

In summary, increasing the diversity or the (constant) forcing term decreases (on 
average) the maximum eigenvalue of the coupled system and thus a Hopf bifurcation 
can be crossed backwards, so that self-oscillations disappear. When applying the 
periodic external forcing on the system formed by self-sustained neurons, coherence 
with the external frequency is difficult to achieve because there is the competing 
effect of mutual neuron synchronization to a different frequency. However, when 
the periodic external forcing is applied on the system of damped neurons, they 
all synchronize to the external forcing, and thus with each other since this is the 
only dynamical regime available to forced damped oscillators (if forcing is not too 
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strong to excite further resonances). Increased coupling strength increases the range 
of uncntraincd self-oscillations. 

Oscillator death by diversity is not particular to this system. In (Mirollo & 
Strogatz 1990) the authors analyze a large system of limit cycle-oscillators with 
mean field coupling and randomly distributed frequencies. They proved that when 
the coupling is sufficiently strong and the distribution of frequencies has a suffi- 
ciently large variance, the system undergoes "amplitude death" . In their approach 
the oscillators pull each other off their limit cycles, which is translated into a stable 
equilibrium point for the coupled system. Thus, this mechanism suggests that the 
quenched noise we introduced in the system "pushes apart" the limit cycles of the 
different neurons, so that their competition enlarges the range of parameters where 
fixed point behavior is stable. 

A qualitative argument explaining the diversity-induced oscillator death in our 
system of coupled neurons goes as follows: We know from Gonze et al. (2005) that 
a single oscillator can switch from a limit cycle to a stable steady state by adding a 
constant mean field (the term containing F in (2.5) but with time-independent F) of 
sufficient strength to Eq. (2.1). A constant light forcing term has the same effect (see 
the zero coupling case in fig. 7). Furthermore we have observed that the amplitude of 
the oscillations decreases with rising diversity (compare figs. 6), but the mean does 
not change. In a system with low diversity we have large oscillations of F around 
that mean value. If this value, taken as a constant, determines a stable steady state, 
then we argue that the large oscillations lead the system into unstable regions, 
whereas, by increasing a the amplitude is decreased and the concentrations do not 
leave neighbourhood of the stable fixed point, thus finding themselves damped all 
the time. This is a possible mechanism for the diversity-induced oscillator death 
phenomenon. 

4. Concluding Remarks 

In this work we have analyzed the role of diversity in favoring the entrainment of a 
system of coupled circadian oscillators. We introduce non-negligible heterogeneity 
in the periods of all neurons in the form of quenched noise. This is achieved by 
rescaling the individual neuronal periods by a scaling factor drawn from a normal 
distribution. The system response to the light-dark cycle periodicity is studied as 
a function of the intcrneuronal coupling strength, external forcing amplitude and 
neuronal heterogeneity. 

Most of the cases of order induced by heterogeneity or noise carried out so far 
(Gammaitoni et al. 1998; Hanggi & Marchesoni 2009; Tessonc et al. 2006, 2007; 
Toral et al. 2009; Pikovsky & Kurths 1997; Ullner et al. 2009), emphasize the fact 
the diversity directly improves oscillator synchronization. In our case the mechanism 
is rather different. Diversity does not improve system synchronization directly. This 
is achieved indirectly, by a leading first to a diversity-induced stabilization of the 
fixed points of the neurons forming the system. Once steady concentrations are 
asymptotically stable, it is much better entrainable by the external forcing, so that 
the damped neurons adapt easily to the external forcing (and then, in addition, they 
appear as synchronized between them). The synchronization arises therefore not as 
a result of a direct diversity-induced neuronal synchronization but indirectly, as a 
result of the diversity-induced oscillator death. Our results indicate therefore that 
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Figure 7. Colour plots of the maximum real part of the average eigenvalues of system of 
Eqs. (2.5-2.9), as a function of a and Lo, at different values of K. Increasing a or increas- 
ing Lo changes this quantity from positive to negative, i.e. transforms the self-sustained 
neurons into damped neurons by stabilizing their constant concentrations fixed points. 
Rising the coupling enlarges the region of self-sustained oscillations. Averaged from 10 
realizations of heterogeneity in 200 neurons. 

the right amount of heterogeneity helps the extended system to respond globally 
in a more coherent way to the external forcing. In addition to the robustness of 
the results against use of non-sinusoidal forcing we have checked that resonances 
in the responses to the external forcing and matching of the circadian period to 
the external forcing appear in more complex models, such as the 10- variable model 
of (Bernard et al. 2007) with diversity in the time scales Tj , or the 4- variable model 
of (Gonze et al. 2005) with heterogeneity in all the reaction rate parameters v^. 
We expect that a similar behavior will be found in models of non-mammalian 
clocks like those of Drosophila (Smolen et al. 2004), Arabidopsis (Locke et al. 2005), 
Neurospora (Hcintzcn & Liu 2007) or Cyanobacteria (Dong & Golden 2008). 

Of course, it is an open question whether the observed diversity in the periods of 
the neurons of the SCN has been tuned by evolution in order to display a maximum 
response to the 24 h dark-light natural cycle. A detailed experimental check of our 
predictions would require to be able to vary the amount of diversity. In this sense 
we suggest the possibility of using chimeric organisms (Low-Zeddies & Takahashi 
2001) to introduce diversity in a controlled way. 
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